#include "function1.h"
#include "TimeIntergratorFactory.h"
#include "VectorFunctionFactory.h"
#include <iomanip>
#include <fstream>
#include <cmath>
#include <ctime>
#include <sys/types.h>
#include <sys/stat.h>

#define TheTimeIntergratorFactory TimeIntergratorFactory<6>::Instance()
#define TheVectorFunctionFactory VectorFunctionFactory<6, 6>::Instance()

int main(int argc, char* argv[])
{
	std::cout << "If the program spends too much time, cut the number of steps of BDF1 in IVP1 by ten percent." << std::endl;
	bool isend = false;
	char c;
	std::string methodId;
	int FunctionId = 1;
	int IVP_Id;
	int p = 0;
	int num_steps;
	int steps;
	double T;
	std::vector<Point<6> > Pnts;
	double x[6];
	std::ifstream data_in("Inputfile");
	TimeIntergrator<6>* Method; 
	VectorFunction<6, 6>* Function;
	int startstep;
	double mid;
	double max_norm[2];
	double E;
	int num_densen;
	double growth1, growth2;
	clock_t beg_time, end_time;

	std::string d = "../results";
	mkdir(d.c_str(), 0777);
	d = "../acctest";
	mkdir(d.c_str(), 0777);
	data_in >> c;
	std::ofstream make_plt("drawfigures.plt");
	make_plt << "unset key" << std::endl;
	make_plt << "set xrange [-1.5: 1.2]" << std::endl;
	while(!isend) {
		if(c == '$')
			std::getline(data_in, methodId, '$');
		data_in >> methodId;
		std::string file_name = "../results/" + methodId;	
		mkdir(file_name.c_str(), 0777);			//create a file to store results.

		if(methodId != "CRK") {					//load and set data.
			data_in >> p;
			Method = TheTimeIntergratorFactory->CreateMethod(methodId, p);
		}
		else {
			Method = TheTimeIntergratorFactory->CreateMethod(methodId);
		}
		data_in >> IVP_Id;
		data_in >> growth1 >> growth2;
		data_in >> num_densen;
		Point<6> P[num_densen];
		int numstep[num_densen];
		double cpu_time[num_densen];
		Function = TheVectorFunctionFactory->CreateFunction(FunctionId);
		data_in >> T >> num_steps;
		for(int i = 0; i != 6; i++) {
			data_in >> x[i];
		}
		
		std::string acclist;
		acclist = "../acctest/acc_" + methodId;
		if(methodId != "CRK")
			acclist = acclist + '_' + std::to_string(p);
		acclist = acclist + "_IVP" + std::to_string(IVP_Id);
		if(c == '&')
			acclist = acclist + "_compair";
		acclist = acclist + ".txt";	//The file name of which stores the error pattern.
		std::ofstream acc_out(acclist);
		acc_out << "num_steps\t";
		if(IVP_Id == 1)
			acc_out << "abs e";
		else
			acc_out << "E_tilde";
		acc_out << "\t\torder\tCPU time" << std::endl;
		
		for(int l = 0; l != num_densen; l++) {
			beg_time = clock();
			Point<6> init_P(x);
			Pnts.push_back(init_P);
			std::string pointlist;
			pointlist = "../results/" + methodId + "/points_";
			if(methodId != "CRK")
				pointlist = pointlist + std::to_string(p) + '_';
			pointlist = pointlist + std::to_string(num_steps) + "_IVP" + std::to_string(IVP_Id);
			if(c == '&')
				pointlist = pointlist + "_compair";
			pointlist = pointlist + ".txt";	//The file name of which stores the points.
			std::ofstream data_out(pointlist);

			//Create the gnuplot script to make plots.
			std::string figures = "../figures";
			mkdir(figures.c_str(), 0777);			//create a file to store figures.
			make_plt << "set title 'Solve IVP" << std::to_string(IVP_Id); 
			make_plt << " with " << std::to_string(num_steps) << "steps by " <<  methodId;
			if(methodId != "CRK")
				make_plt << std::to_string(p);
			make_plt << "'" << std::endl;
			make_plt << "set term post color" << std::endl;
			make_plt << "set output '../figures/" << methodId << "_" << std::to_string(num_steps);
			if(methodId != "CRK")
				make_plt << "_" << std::to_string(p);
			make_plt << "_IVP" << std::to_string(IVP_Id);
			if(c == '&')
				make_plt << "_compair";
			make_plt << ".eps'" << std::endl;
			make_plt << "plot \"" << pointlist << "\" using 1:2 w l lt 3 lw 1" << std::endl;

			for(int i = 0; i != 2; i++) {
				data_out << std::fixed << std::setprecision(15) << init_P[i] << '\t';
			}
			data_out << std::endl;
			startstep = 0;
			steps = 0;
			if(methodId != "CRK") {
				Method->one_step_method(Pnts, T/num_steps, *Function);
				for(auto i = Pnts.begin(); i != Pnts.end(); i++) {
					startstep++;
					steps++;
					for(int j = 0; j != 2; j++) {
						data_out << (*i)[j] << '\t';
					}
					data_out << std::endl;
				}
			}
			std::cout << "start:" << methodId;
			if(methodId != "CRK")
				std::cout << "_" << std::to_string(p);
			std::cout << "_IVP" << std::to_string(IVP_Id) << std::endl;
			for(int i = 0; i != num_steps - startstep; i++) {
				Method->one_step_method(Pnts, T/num_steps, *Function);
				steps++;
				for(int j = 0; j != 2; j++) {
					data_out << (*(Pnts.end() - 1))[j] << '\t';
				}
				data_out << std::endl;
			}
			end_time = clock();
			cpu_time[l] = double(end_time - beg_time)/CLOCKS_PER_SEC;
			std::cout << *(Pnts.end() - 1) << std::endl;
			std::cout << steps << std::endl;
			numstep[l] = num_steps;
			num_steps = num_steps*growth1/growth2;
			data_out.close();
			P[l] = *(Pnts.end() - 1);		//store the end points.
			Pnts.clear();
		}

		//Accuracy and convergence rates tests.
		for(int l = 0; l != num_densen - 1; l++) {
			max_norm[1] = 0;
			for(int i = 0; i != 6; i++) {
				if(IVP_Id == 2)
					mid = P[l + 1][i] - P[l][i];
				else
					mid = x[i] - P[l][i];
				if(max_norm[1] < abs(mid)) 
					max_norm[1] = abs(mid);
			}
			acc_out << numstep[l + 1] << "\t" << max_norm[1] << '\t';
			if(l != 0) {
				E = log2(max_norm[0]/max_norm[1])/log2((double)growth1/growth2);
				acc_out << E << '\t' << cpu_time[l] << std::endl;
			}
			else {
				acc_out << std::endl;
			}
			max_norm[0] = max_norm[1];
		}
		acc_out.close();
		data_in >> c;
		if(c == '#') {
			isend = true;
		}
	}

	make_plt << "shell";
	make_plt.close();
	data_in.close();
}

